library(spacc)
library(data.table)
library(ggplot2)
data_path <- "../data/"
header <- fread(file.path(data_path, "austria_header.csv"))
species <- fread(file.path(data_path, "austria_species.csv"))
cat(sprintf("Plots: %d\nSpecies records: %d\nUnique species: %d\n",
nrow(header), nrow(species), uniqueN(species$WFO_TAXON)))
## Plots: 52526
## Species records: 852075
## Unique species: 2700
species[, .N, by = STATUS]
## STATUS N
## <char> <int>
## 1: native 826055
## 2: neo 26020
Split by native/alien status and create separate site-by-species matrices sharing the same row order. We also assign each plot to a decade for temporal analysis.
rds_mats <- file.path(comp_dir, "matrices.rds")
# Keep only plots present in both header and species
common_ids <- intersect(header$PlotObservationID, species$PlotObservationID)
header <- header[PlotObservationID %in% common_ids]
# Assign decade
header[, decade := ifelse(Year < 1970, "pre-1970",
paste0(floor(Year / 10) * 10, "s"))]
# Coordinates
coords <- header[, .(x = Longitude, y = Latitude)]
# Helper: long-form species table -> wide presence/absence matrix
make_pa_matrix <- function(sp_data, plot_ids) {
sp_data <- sp_data[PlotObservationID %in% plot_ids]
wide <- dcast(sp_data, PlotObservationID ~ WFO_TAXON,
fun.aggregate = function(x) as.integer(length(x) > 0),
value.var = "WFO_TAXON")
wide <- wide[match(plot_ids, PlotObservationID)]
mat <- as.matrix(wide[, -1])
mat[is.na(mat)] <- 0L
rownames(mat) <- wide$PlotObservationID
mat
}
if (file.exists(rds_mats)) {
cat("Loading cached matrices...\n")
mats <- readRDS(rds_mats)
native_mat <- mats$native; alien_mat <- mats$alien; all_mat <- mats$all
rm(mats)
} else {
plot_ids <- header$PlotObservationID
native_mat <- make_pa_matrix(species[STATUS == "native"], plot_ids)
alien_mat <- make_pa_matrix(species[STATUS == "neo"], plot_ids)
all_mat <- make_pa_matrix(species, plot_ids)
saveRDS(list(native = native_mat, alien = alien_mat, all = all_mat), rds_mats)
}
## Loading cached matrices...
cat(sprintf("Native species: %d\nAlien species: %d\nAll species: %d\n",
ncol(native_mat), ncol(alien_mat), ncol(all_mat)))
## Native species: 2378
## Alien species: 323
## All species: 2700
# Decade summary
cat("\nPlots per decade:\n")
##
## Plots per decade:
print(header[, .N, by = decade][order(decade)])
## decade N
## <char> <int>
## 1: 1970s 2386
## 2: 1980s 3171
## 3: 1990s 8814
## 4: 2000s 4992
## 5: 2010s 8285
## 6: 2020s 104
## 7: pre-1970 2105
By splitting into decades we get smaller matrices (1K-9K sites each), enabling high-seed-count runs with temporal resolution.
decades <- sort(unique(header$decade))
rds_decade <- file.path(comp_dir, "decade_results.rds")
if (file.exists(rds_decade)) {
cat("Loading cached decade results...\n")
decade_results <- readRDS(rds_decade)
} else {
decade_results <- lapply(decades, function(d) {
idx <- which(header$decade == d)
n <- length(idx)
cat(sprintf(" %s: %d sites\n", d, n))
spacc(all_mat[idx, , drop = FALSE], coords[idx, ],
n_seeds = min(100, n), method = "knn",
distance = "haversine", seed = 42)
})
names(decade_results) <- decades
saveRDS(decade_results, rds_decade)
}
## Loading cached decade results...
decade_combined <- do.call(c, decade_results)
plot(decade_combined) + ggtitle("Spatial accumulation by decade")
Same spatial walks, different species tallies. Using the
groups argument for efficiency.
status_vec <- species[match(colnames(all_mat), WFO_TAXON), STATUS]
rds_decade_na <- file.path(comp_dir, "decade_na_results.rds")
if (file.exists(rds_decade_na)) {
cat("Loading cached decade native/alien results...\n")
decade_na_results <- readRDS(rds_decade_na)
} else {
decade_na_results <- lapply(decades, function(d) {
idx <- which(header$decade == d)
n <- length(idx)
spacc(all_mat[idx, , drop = FALSE], coords[idx, ],
n_seeds = min(100, n), method = "knn",
distance = "haversine", groups = status_vec, seed = 42)
})
names(decade_na_results) <- decades
saveRDS(decade_na_results, rds_decade_na)
}
## Loading cached decade native/alien results...
for (d in decades) {
print(plot_grouped_normalized(decade_na_results[[d]],
title = sprintf("%s: Native vs Alien (normalized)", d)))
}
For overall comparisons we use the full dataset with moderate seed counts.
rds_pooled <- file.path(comp_dir, "pooled.rds")
if (file.exists(rds_pooled)) {
cat("Loading cached pooled results...\n")
p <- readRDS(rds_pooled)
sac_native <- p$native; sac_alien <- p$alien; sac_all <- p$all
rm(p)
} else {
sac_native <- spacc(native_mat, coords, n_seeds = 30, method = "knn",
distance = "haversine", seed = 42)
sac_alien <- spacc(alien_mat, coords, n_seeds = 30, method = "knn",
distance = "haversine", seed = 42)
sac_all <- spacc(all_mat, coords, n_seeds = 30, method = "knn",
distance = "haversine", seed = 42)
saveRDS(list(native = sac_native, alien = sac_alien, all = sac_all), rds_pooled)
}
## Loading cached pooled results...
plot_normalized(Native = sac_native, Alien = sac_alien,
title = "Pooled: Native vs Alien (normalized)")
rds_compare <- file.path(comp_dir, "compare.rds")
if (file.exists(rds_compare)) {
comp <- readRDS(rds_compare)
} else {
comp <- compare(sac_native, sac_alien, method = "permutation", n_perm = 999)
saveRDS(comp, rds_compare)
}
print(comp)
## Comparison: sac_native vs sac_alien
## ----------------------------------------
## Method: permutation (n=999)
## AUC difference: 50036933.1 (p = 0.000***)
## Saturation: sac_native at 18661 sites, sac_alien at 24060 sites
##
## sac_native saturates faster.
plot(comp)
Map per-site species richness across Austria to reveal spatial hotspots. These are the real data equivalent of the simulated heatmaps from Day 2 theory.
library(sf)
library(rnaturalearth)
austria_sf <- ne_countries(scale = "medium", country = "Austria", returnclass = "sf")
# Per-site richness from the matrices
site_richness <- data.frame(
x = coords$x,
y = coords$y,
native = rowSums(native_mat),
alien = rowSums(alien_mat),
total = rowSums(all_mat)
)
site_richness$alien_pct <- site_richness$alien / site_richness$total * 100
# Native heatmap
p_native <- ggplot() +
geom_sf(data = austria_sf, fill = "grey95", color = "#666666", linewidth = 0.8) +
geom_point(data = site_richness, aes(x = x, y = y, color = native),
size = 0.4, alpha = 0.6) +
geom_sf(data = austria_sf, fill = NA, color = "#333333", linewidth = 1) +
scale_color_gradient(low = "#E8F5E9", high = "#1B5E20", name = "Native\nSpecies") +
labs(title = "Native Species Richness per Plot",
x = "Longitude", y = "Latitude") +
theme_minimal(base_size = 12) +
coord_sf(xlim = c(9.3, 17.2), ylim = c(46.3, 49.1))
# Alien heatmap
p_alien <- ggplot() +
geom_sf(data = austria_sf, fill = "grey95", color = "#666666", linewidth = 0.8) +
geom_point(data = site_richness, aes(x = x, y = y, color = alien),
size = 0.4, alpha = 0.6) +
geom_sf(data = austria_sf, fill = NA, color = "#333333", linewidth = 1) +
scale_color_gradient(low = "#FFEBEE", high = "#B71C1C", name = "Alien\nSpecies") +
labs(title = "Alien Species Richness per Plot",
x = "Longitude", y = "Latitude") +
theme_minimal(base_size = 12) +
coord_sf(xlim = c(9.3, 17.2), ylim = c(46.3, 49.1))
gridExtra::grid.arrange(p_native, p_alien, ncol = 1)
# Proportion alien map
ggplot() +
geom_sf(data = austria_sf, fill = "grey95", color = "#666666", linewidth = 0.8) +
geom_point(data = site_richness, aes(x = x, y = y, color = alien_pct),
size = 0.4, alpha = 0.6) +
geom_sf(data = austria_sf, fill = NA, color = "#333333", linewidth = 1) +
scale_color_gradient2(low = "#1B5E20", mid = "#FFF9C4", high = "#B71C1C",
midpoint = median(site_richness$alien_pct, na.rm = TRUE),
name = "Alien %") +
labs(title = "Alien Species as % of Total Richness",
x = "Longitude", y = "Latitude") +
theme_minimal(base_size = 12) +
coord_sf(xlim = c(9.3, 17.2), ylim = c(46.3, 49.1))
Accumulation curves reveal where species are concentrated. By starting spatial walks from different regions, we can detect invasion hotspots: if aliens accumulate faster when starting near cities, that confirms urban areas are introduction points. Conversely, native species should accumulate faster from alpine starts. We pick seeds from three regions: the Vienna lowlands (east, urban), the Alps (west, natural), and central Austria.
# Define three geographic regions
regions <- list(
"Vienna & East" = list(lon = c(15.5, 17.0), lat = c(47.5, 48.5), color = "#D32F2F"),
"Alps (West)" = list(lon = c(10.0, 12.5), lat = c(46.5, 47.5), color = "#1976D2"),
"Center" = list(lon = c(13.0, 15.0), lat = c(47.5, 48.5), color = "#7B1FA2")
)
# Find sites in each region
region_indices <- lapply(regions, function(r) {
which(coords$x >= r$lon[1] & coords$x <= r$lon[2] &
coords$y >= r$lat[1] & coords$y <= r$lat[2])
})
# Map showing regions and plot locations
region_df <- do.call(rbind, lapply(names(regions), function(nm) {
idx <- region_indices[[nm]]
data.frame(x = coords$x[idx], y = coords$y[idx], region = nm)
}))
cat("Sites per region:\n")
## Sites per region:
for (nm in names(regions)) cat(sprintf(" %s: %d sites\n", nm, length(region_indices[[nm]])))
## Vienna & East: 10810 sites
## Alps (West): 1732 sites
## Center: 5871 sites
ggplot() +
geom_sf(data = austria_sf, fill = "grey95", color = "#666666", linewidth = 0.8) +
geom_point(data = data.frame(x = coords$x, y = coords$y),
aes(x = x, y = y), color = "grey80", size = 0.2, alpha = 0.3) +
geom_point(data = region_df, aes(x = x, y = y, color = region), size = 0.8, alpha = 0.6) +
geom_sf(data = austria_sf, fill = NA, color = "#333333", linewidth = 1) +
scale_color_manual(values = c("Vienna & East" = "#D32F2F",
"Alps (West)" = "#1976D2",
"Center" = "#7B1FA2")) +
labs(title = "Starting Regions for Invasion Hotspot Detection",
subtitle = "Where you start reveals where species are concentrated",
x = "Longitude", y = "Latitude", color = "Region") +
theme_minimal(base_size = 12) +
coord_sf(xlim = c(9.3, 17.2), ylim = c(46.3, 49.1))
rds_startpoint <- file.path(comp_dir, "startpoint.rds")
if (file.exists(rds_startpoint)) {
sp_results <- readRDS(rds_startpoint)
} else {
n_seeds_per_region <- 5
set.seed(42)
nc <- parallel::detectCores(logical = FALSE)
sp_results <- lapply(names(regions), function(nm) {
idx <- region_indices[[nm]]
# Sample seed sites (0-based for C++)
seeds <- sample(idx, min(n_seeds_per_region, length(idx)), replace = FALSE) - 1L
# Run from these specific seeds on ALL sites (full dataset kNN walk)
curves_all <- spacc:::cpp_knn_kdtree_parallel_seeds(
all_mat, coords$x, coords$y, seeds, n_cores = nc,
progress = FALSE, distance = "haversine")
curves_native <- spacc:::cpp_knn_kdtree_parallel_seeds(
native_mat, coords$x, coords$y, seeds, n_cores = nc,
progress = FALSE, distance = "haversine")
curves_alien <- spacc:::cpp_knn_kdtree_parallel_seeds(
alien_mat, coords$x, coords$y, seeds, n_cores = nc,
progress = FALSE, distance = "haversine")
list(all = curves_all, native = curves_native, alien = curves_alien,
seeds = seeds, n_seeds = length(seeds))
})
names(sp_results) <- names(regions)
saveRDS(sp_results, rds_startpoint)
}
# Plot: all-species curves by starting region
curves_df <- do.call(rbind, lapply(names(sp_results), function(nm) {
mat <- sp_results[[nm]]$all
n <- ncol(mat)
total <- max(mat)
data.frame(
sites = seq_len(n),
mean = colMeans(mat) / total,
region = nm
)
}))
curves_df$region <- factor(curves_df$region, levels = names(regions))
p1 <- ggplot(curves_df, aes(x = sites, y = mean, color = region)) +
geom_line(linewidth = 1) +
scale_color_manual(values = c("Vienna & East" = "#D32F2F",
"Alps (West)" = "#1976D2",
"Center" = "#7B1FA2")) +
scale_y_continuous(labels = scales::percent) +
labs(title = "All Species: Curves by Starting Region",
x = "Sites sampled", y = "Proportion of total species",
color = "Start Region") +
theme_minimal(base_size = 12)
# Plot: native curves by starting region
native_df <- do.call(rbind, lapply(names(sp_results), function(nm) {
mat <- sp_results[[nm]]$native
n <- ncol(mat)
total <- max(mat)
data.frame(
sites = seq_len(n),
mean = colMeans(mat) / total,
region = nm
)
}))
native_df$region <- factor(native_df$region, levels = names(regions))
p2 <- ggplot(native_df, aes(x = sites, y = mean, color = region)) +
geom_line(linewidth = 1) +
scale_color_manual(values = c("Vienna & East" = "#D32F2F",
"Alps (West)" = "#1976D2",
"Center" = "#7B1FA2")) +
scale_y_continuous(labels = scales::percent) +
labs(title = "Native Species: Curves by Starting Region",
x = "Sites sampled", y = "Proportion of total native species",
color = "Start Region") +
theme_minimal(base_size = 12)
# Plot: alien curves by starting region
alien_df <- do.call(rbind, lapply(names(sp_results), function(nm) {
mat <- sp_results[[nm]]$alien
n <- ncol(mat)
total <- max(mat)
data.frame(
sites = seq_len(n),
mean = colMeans(mat) / total,
region = nm
)
}))
alien_df$region <- factor(alien_df$region, levels = names(regions))
p3 <- ggplot(alien_df, aes(x = sites, y = mean, color = region)) +
geom_line(linewidth = 1) +
scale_color_manual(values = c("Vienna & East" = "#D32F2F",
"Alps (West)" = "#1976D2",
"Center" = "#7B1FA2")) +
scale_y_continuous(labels = scales::percent) +
labs(title = "Alien Species: Curves by Starting Region",
x = "Sites sampled", y = "Proportion of total alien species",
color = "Start Region") +
theme_minimal(base_size = 12)
gridExtra::grid.arrange(p1, p2, p3, ncol = 1)
# Zoom into first 1000 sites to see early differences clearly
n_zoom <- min(1000, ncol(sp_results[[1]]$all))
early_df <- do.call(rbind, lapply(names(sp_results), function(nm) {
mat_n <- sp_results[[nm]]$native
mat_a <- sp_results[[nm]]$alien
total_n <- max(mat_n)
total_a <- max(mat_a)
data.frame(
sites = rep(seq_len(n_zoom), 2),
mean = c(colMeans(mat_n)[1:n_zoom] / total_n,
colMeans(mat_a)[1:n_zoom] / total_a),
group = rep(c("Native", "Alien"), each = n_zoom),
region = nm
)
}))
early_df$label <- paste(early_df$region, "-", early_df$group)
early_df$region <- factor(early_df$region, levels = names(regions))
ggplot(early_df, aes(x = sites, y = mean, color = region, linetype = group)) +
geom_line(linewidth = 0.8) +
scale_color_manual(values = c("Vienna & East" = "#D32F2F",
"Alps (West)" = "#1976D2",
"Center" = "#7B1FA2")) +
scale_linetype_manual(values = c("Native" = "solid", "Alien" = "dashed")) +
scale_y_continuous(labels = scales::percent) +
labs(title = sprintf("Early Accumulation (first %d sites): Native vs Alien by Region", n_zoom),
subtitle = "Faster alien accumulation from Vienna confirms urban invasion hotspots",
x = "Sites sampled", y = "Proportion of total species",
color = "Start Region", linetype = "Species Group") +
theme_minimal(base_size = 12)
Fit asymptotic models to estimate total richness for each group.
rds_extrap <- file.path(comp_dir, "extrapolate.rds")
if (file.exists(rds_extrap)) {
e <- readRDS(rds_extrap)
fit_native <- e$native; fit_alien <- e$alien
rm(e)
} else {
fit_native <- extrapolate(sac_native, model = "michaelis-menten")
fit_alien <- extrapolate(sac_alien, model = "michaelis-menten")
saveRDS(list(native = fit_native, alien = fit_alien), rds_extrap)
}
cat("Native:\n"); print(fit_native)
## Native:
## Extrapolation: michaelis-menten
## ------------------------------
## Estimated asymptote: 2525.4 species
## 95% CI: 2522.4 - 2528.5
## AIC: 354656.5
## Observed: 2378.0 species (94% of estimated)
cat("\nAlien:\n"); print(fit_alien)
##
## Alien:
## Extrapolation: michaelis-menten
## ------------------------------
## Estimated asymptote: 410.9 species
## 95% CI: 410.0 - 411.8
## AIC: 233579.6
## Observed: 323.0 species (79% of estimated)
plot(fit_native)
plot(fit_alien)
For each hexagon in an equal-area grid, run a spatial accumulation with 10 seeds on the within-hex plots, then fit Michaelis-Menten to estimate total species richness (S_max). This gives estimated richness per hexagon rather than raw observed counts.
library(hexify)
library(sf)
library(rnaturalearth)
austria_sf <- ne_countries(scale = "medium", country = "Austria", returnclass = "sf")
grid <- hex_grid(area_km2 = 750)
# Assign plots to hexagons using hexify()
hx <- hexify(data.frame(lon = header$Longitude, lat = header$Latitude), grid = grid)
plot_hex <- data.table(
PlotObservationID = header$PlotObservationID,
lon = header$Longitude,
lat = header$Latitude,
cell_id = hx[["cell_id"]]
)
# Count plots per hex
hex_counts <- plot_hex[, .N, by = cell_id]
# Keep hexagons with >= 20 plots (need enough for meaningful SAC + extrapolation)
keep_cells <- hex_counts[N >= 20]$cell_id
cat(sprintf("Hexagons with >= 20 plots: %d (out of %d)\n",
length(keep_cells), nrow(hex_counts)))
## Hexagons with >= 20 plots: 118 (out of 130)
rds_hex_sac <- file.path(comp_dir, "hex_sac_estimates.rds")
if (file.exists(rds_hex_sac)) {
cat("Loading cached per-hexagon SAC estimates...\n")
hex_estimates <- readRDS(rds_hex_sac)
} else {
hex_estimates <- rbindlist(lapply(keep_cells, function(cid) {
idx <- which(plot_hex$cell_id == cid)
n <- length(idx)
# Run spacc within this hexagon
result <- tryCatch({
sac_nat <- spacc(native_mat[idx, , drop = FALSE],
coords[idx, ], n_seeds = 10, method = "knn",
distance = "haversine", seed = 42)
sac_ali <- spacc(alien_mat[idx, , drop = FALSE],
coords[idx, ], n_seeds = 10, method = "knn",
distance = "haversine", seed = 42)
sac_tot <- spacc(all_mat[idx, , drop = FALSE],
coords[idx, ], n_seeds = 10, method = "knn",
distance = "haversine", seed = 42)
# Extrapolate Michaelis-Menten
fit_nat <- tryCatch(extrapolate(sac_nat, model = "michaelis-menten"),
error = function(e) NULL)
fit_ali <- tryCatch(extrapolate(sac_ali, model = "michaelis-menten"),
error = function(e) NULL)
fit_tot <- tryCatch(extrapolate(sac_tot, model = "michaelis-menten"),
error = function(e) NULL)
# Observed species counts
obs_native <- ncol(native_mat[idx, colSums(native_mat[idx, , drop = FALSE]) > 0,
drop = FALSE])
obs_alien <- ncol(alien_mat[idx, colSums(alien_mat[idx, , drop = FALSE]) > 0,
drop = FALSE])
obs_total <- ncol(all_mat[idx, colSums(all_mat[idx, , drop = FALSE]) > 0,
drop = FALSE])
data.table(
cell_id = cid,
n_plots = n,
obs_native = obs_native,
obs_alien = obs_alien,
obs_total = obs_total,
est_native = if (!is.null(fit_nat)) fit_nat$asymptote else NA_real_,
est_alien = if (!is.null(fit_ali)) fit_ali$asymptote else NA_real_,
est_total = if (!is.null(fit_tot)) fit_tot$asymptote else NA_real_
)
}, error = function(e) {
data.table(cell_id = cid, n_plots = n,
obs_native = NA, obs_alien = NA, obs_total = NA,
est_native = NA, est_alien = NA, est_total = NA)
})
cat(sprintf(" hex %s: %d plots, est_native=%.0f, est_alien=%.0f\n",
cid, n,
result$est_native, result$est_alien))
result
}))
saveRDS(hex_estimates, rds_hex_sac)
}
## Loading cached per-hexagon SAC estimates...
cat(sprintf("Estimated richness for %d hexagons\n", nrow(hex_estimates)))
## Estimated richness for 118 hexagons
cat(sprintf(" Native range: %.0f - %.0f\n",
min(hex_estimates$est_native, na.rm = TRUE),
max(hex_estimates$est_native, na.rm = TRUE)))
## Native range: 270 - 3614
cat(sprintf(" Alien range: %.0f - %.0f\n",
min(hex_estimates$est_alien, na.rm = TRUE),
max(hex_estimates$est_alien, na.rm = TRUE)))
## Alien range: 2 - 872
hex_estimates[, completeness_native := obs_native / est_native * 100]
hex_estimates[, completeness_alien := obs_alien / est_alien * 100]
hex_estimates[, completeness_total := obs_total / est_total * 100]
cat(sprintf("Median completeness: native=%.1f%%, alien=%.1f%%\n",
median(hex_estimates$completeness_native, na.rm = TRUE),
median(hex_estimates$completeness_alien, na.rm = TRUE)))
## Median completeness: native=64.8%, alien=56.8%
hex_polys <- cell_to_sf(hex_estimates$cell_id, grid = grid)
hex_polys <- merge(hex_polys, hex_estimates, by = "cell_id")
# Estimated native richness
p1 <- ggplot() +
geom_sf(data = austria_sf, fill = "grey95", color = "#333333", linewidth = 0.8) +
geom_sf(data = hex_polys, aes(fill = est_native), color = "white", linewidth = 0.3) +
geom_sf(data = austria_sf, fill = NA, color = "#333333", linewidth = 1) +
scale_fill_gradient(low = "white", high = "#2E7D32",
name = "Estimated\nnative S_max") +
labs(title = "Estimated Native Species Richness (Michaelis-Menten S_max per hexagon)",
x = "Longitude", y = "Latitude") +
theme_minimal(base_size = 12) +
coord_sf(xlim = c(9.3, 17.2), ylim = c(46.3, 49.1))
# Estimated alien richness
p2 <- ggplot() +
geom_sf(data = austria_sf, fill = "grey95", color = "#333333", linewidth = 0.8) +
geom_sf(data = hex_polys, aes(fill = est_alien), color = "white", linewidth = 0.3) +
geom_sf(data = austria_sf, fill = NA, color = "#333333", linewidth = 1) +
scale_fill_gradient(low = "white", high = "#D32F2F",
name = "Estimated\nalien S_max") +
labs(title = "Estimated Alien Species Richness (Michaelis-Menten S_max per hexagon)",
x = "Longitude", y = "Latitude") +
theme_minimal(base_size = 12) +
coord_sf(xlim = c(9.3, 17.2), ylim = c(46.3, 49.1))
# Completeness: % observed / estimated
p3 <- ggplot() +
geom_sf(data = austria_sf, fill = "grey95", color = "#333333", linewidth = 0.8) +
geom_sf(data = hex_polys, aes(fill = completeness_alien), color = "white", linewidth = 0.3) +
geom_sf(data = austria_sf, fill = NA, color = "#333333", linewidth = 1) +
scale_fill_distiller(palette = "YlOrRd", direction = 1,
name = "Alien\ncompleteness %") +
labs(title = "Alien Sampling Completeness per Hexagon (observed / estimated)",
x = "Longitude", y = "Latitude") +
theme_minimal(base_size = 12) +
coord_sf(xlim = c(9.3, 17.2), ylim = c(46.3, 49.1))
gridExtra::grid.arrange(p1, p2, p3, ncol = 1)
Same approach but grouping by EUNIS habitat type rather than hexagon.
# Define habitat variables (also used later in section 16)
habitat <- header$Eunis_lvl1
hab_counts <- table(habitat)
keep_habs <- names(hab_counts[hab_counts >= 100])
keep_habs <- keep_habs[!keep_habs %in% c("", "~")]
rds_hab_sac <- file.path(comp_dir, "habitat_sac_estimates.rds")
if (file.exists(rds_hab_sac)) {
cat("Loading cached per-habitat SAC estimates...\n")
hab_estimates <- readRDS(rds_hab_sac)
} else {
hab_estimates <- rbindlist(lapply(keep_habs, function(h) {
idx <- which(habitat == h)
n <- length(idx)
result <- tryCatch({
sac_nat <- spacc(native_mat[idx, , drop = FALSE],
coords[idx, ], n_seeds = 10, method = "knn",
distance = "haversine", seed = 42)
sac_ali <- spacc(alien_mat[idx, , drop = FALSE],
coords[idx, ], n_seeds = 10, method = "knn",
distance = "haversine", seed = 42)
fit_nat <- tryCatch(extrapolate(sac_nat, model = "michaelis-menten"),
error = function(e) NULL)
fit_ali <- tryCatch(extrapolate(sac_ali, model = "michaelis-menten"),
error = function(e) NULL)
obs_native <- ncol(native_mat[idx, colSums(native_mat[idx, , drop = FALSE]) > 0,
drop = FALSE])
obs_alien <- ncol(alien_mat[idx, colSums(alien_mat[idx, , drop = FALSE]) > 0,
drop = FALSE])
data.table(
habitat = h,
n_plots = n,
obs_native = obs_native,
obs_alien = obs_alien,
est_native = if (!is.null(fit_nat)) fit_nat$asymptote else NA_real_,
est_alien = if (!is.null(fit_ali)) fit_ali$asymptote else NA_real_,
completeness_native = if (!is.null(fit_nat)) obs_native / fit_nat$asymptote * 100 else NA,
completeness_alien = if (!is.null(fit_ali)) obs_alien / fit_ali$asymptote * 100 else NA
)
}, error = function(e) {
data.table(habitat = h, n_plots = n,
obs_native = NA, obs_alien = NA,
est_native = NA, est_alien = NA,
completeness_native = NA, completeness_alien = NA)
})
cat(sprintf(" %s: %d plots, est_native=%.0f, est_alien=%.0f\n",
h, n, result$est_native, result$est_alien))
result
}))
saveRDS(hab_estimates, rds_hab_sac)
}
## Loading cached per-habitat SAC estimates...
print(hab_estimates)
## habitat n_plots obs_native obs_alien est_native est_alien
## <char> <int> <int> <int> <num> <num>
## 1: P 223 287 2 430.150 4.722158
## 2: Q 1102 761 43 1059.577 87.765180
## 3: R 14193 2100 174 2242.250 263.673285
## 4: S 1109 1361 67 1635.407 1099.594735
## 5: T 7001 1583 135 1803.650 201.062973
## 6: U 327 834 56 1456.727 560.110861
## 7: V 4642 875 195 1096.790 280.025930
## completeness_native completeness_alien
## <num> <num>
## 1: 66.72092 42.353520
## 2: 71.82113 48.994373
## 3: 93.65592 65.990758
## 4: 83.22087 6.093154
## 5: 87.76647 67.143143
## 6: 57.25165 9.998021
## 7: 79.77827 69.636408
Compare kNN (spatially structured) vs random (null model) accumulation. We use the 1990s decade (largest chunk) for speed.
idx_90s <- which(header$decade == "1990s")
rds_methods <- file.path(comp_dir, "methods_comparison.rds")
if (file.exists(rds_methods)) {
m <- readRDS(rds_methods)
sac_knn_m <- m$knn; sac_random_m <- m$random
rm(m)
} else {
sac_knn_m <- spacc(all_mat[idx_90s, , drop = FALSE], coords[idx_90s, ],
method = "knn", n_seeds = 50,
distance = "haversine", seed = 42)
sac_random_m <- spacc(all_mat[idx_90s, , drop = FALSE], coords[idx_90s, ],
method = "random", n_seeds = 50,
distance = "haversine", seed = 42)
saveRDS(list(knn = sac_knn_m, random = sac_random_m), rds_methods)
}
methods_combined <- c(kNN = sac_knn_m, Random = sac_random_m)
plot(methods_combined) + ggtitle("kNN vs Random (1990s subset)")
Track effective species diversity (q = 0, 1, 2) as sites accumulate. Using 1990s decade.
rds_hill <- file.path(comp_dir, "hill.rds")
if (file.exists(rds_hill)) {
hill <- readRDS(rds_hill)
} else {
hill <- spaccHill(all_mat[idx_90s, , drop = FALSE], coords[idx_90s, ],
q = c(0, 1, 2), n_seeds = 30,
method = "knn", distance = "haversine", seed = 42)
saveRDS(hill, rds_hill)
}
print(hill)
## spacc Hill numbers: 8814 sites, 2700 species, 30 seeds
## Orders (q): 0, 1, 2
## Method: knn
plot(hill)
Compare Hill number accumulation separately for native and alien species.
rds_hill_na <- file.path(comp_dir, "hill_native_alien.rds")
if (file.exists(rds_hill_na)) {
hill_na <- readRDS(rds_hill_na)
} else {
hill_native <- spaccHill(native_mat[idx_90s, , drop = FALSE], coords[idx_90s, ],
q = c(0, 1, 2), n_seeds = 30,
method = "knn", distance = "haversine", seed = 42)
hill_alien <- spaccHill(alien_mat[idx_90s, , drop = FALSE], coords[idx_90s, ],
q = c(0, 1, 2), n_seeds = 30,
method = "knn", distance = "haversine", seed = 42)
hill_na <- list(native = hill_native, alien = hill_alien)
saveRDS(hill_na, rds_hill_na)
}
cat("Native Hill:", class(hill_na$native), "\n")
## Native Hill: spacc_hill
cat("Alien Hill:", class(hill_na$alien), "\n")
## Alien Hill: spacc_hill
Track how turnover and nestedness change as the sampling window expands.
rds_beta <- file.path(comp_dir, "beta.rds")
if (file.exists(rds_beta)) {
beta <- readRDS(rds_beta)
} else {
beta <- spaccBeta(all_mat[idx_90s, , drop = FALSE], coords[idx_90s, ],
n_seeds = 30, index = "sorensen",
distance = "haversine", seed = 42)
saveRDS(beta, rds_beta)
}
print(beta)
## spacc beta diversity: 8814 sites, 30 seeds
## Index: sorensen, Method: knn
## Mean final beta: 0.976 (turnover: 0.000, nestedness: 0.976)
plot(beta)
How many sites do we need for 90%, 95%, 99% completeness?
rds_cov <- file.path(comp_dir, "coverage.rds")
if (file.exists(rds_cov)) {
cov <- readRDS(rds_cov)
} else {
cov <- spaccCoverage(all_mat[idx_90s, , drop = FALSE], coords[idx_90s, ],
n_seeds = 30,
distance = "haversine", seed = 42)
saveRDS(cov, rds_cov)
}
print(cov)
## spacc coverage: 8814 sites, 2700 species, 30 seeds
## Mean final coverage: 99.9%
## Mean final richness: 2053.0 species
plot(cov)
# Interpolate to standard coverage targets
interp <- interpolateCoverage(cov, target = c(0.90, 0.95, 0.99))
cat("\nSites needed for target coverage:\n")
##
## Sites needed for target coverage:
print(round(colMeans(interp, na.rm = TRUE), 1))
## C90 C95 C99
## 134.1 192.9 463.1
Compute accumulation metrics for every site as a starting point. Using 1990s subset.
rds_metrics <- file.path(comp_dir, "metrics.rds")
if (file.exists(rds_metrics)) {
metrics <- readRDS(rds_metrics)
} else {
metrics <- spaccMetrics(all_mat[idx_90s, , drop = FALSE], coords[idx_90s, ],
metrics = c("slope_10", "half_richness", "auc"),
method = "knn", distance = "haversine")
saveRDS(metrics, rds_metrics)
}
print(metrics)
## spacc_metrics: 8814 sites, 2700 species
## Method: knn
## Metrics: slope_10, half_richness, auc
summary(metrics)
## Metric summary:
## slope_10: mean=0.82, sd=0.17, range=[0.47, 1.23]
## half_richness: mean=2050.14, sd=486.35, range=[1204.00, 3454.00]
## auc: mean=1585.70, sd=38.80, range=[1503.78, 1658.52]
plot(metrics, type = "points")
Decompose regional diversity into alpha, beta, and gamma components.
rds_part <- file.path(comp_dir, "partition.rds")
if (file.exists(rds_part)) {
partition <- readRDS(rds_part)
} else {
partition <- diversityPartition(all_mat[idx_90s, , drop = FALSE], q = c(0, 1, 2),
coords = data.frame(x = coords$x[idx_90s], y = coords$y[idx_90s]))
saveRDS(partition, rds_part)
}
print(partition)
## Alpha-Beta-Gamma Diversity Partitioning
## 8814 sites, 2700 species
##
## q alpha beta gamma
## 0 30.54 67.22 2053.00
## 1 26.54 29.17 774.29
## 2 21.09 23.27 490.89
##
## Interpretation:
## Alpha = mean effective species per site
## Beta = effective number of communities (1 to n_sites)
## Gamma = regional effective species (gamma = alpha x beta)
plot(partition)
rds_alpha <- file.path(comp_dir, "alpha.rds")
if (file.exists(rds_alpha)) {
alpha <- readRDS(rds_alpha)
} else {
alpha <- alphaDiversity(all_mat[idx_90s, , drop = FALSE], q = c(0, 1, 2),
coords = data.frame(x = coords$x[idx_90s], y = coords$y[idx_90s]))
saveRDS(alpha, rds_alpha)
}
print(alpha)
## spacc alpha diversity: 8814 sites, 2700 species
## Orders (q): 0, 1, 2
## Coordinates: available
plot(alpha)
Compare simulation-based curves to exact analytical expectations. Using 1990s subset.
cole <- coleman(all_mat[idx_90s, , drop = FALSE])
mt <- mao_tau(all_mat[idx_90s, , drop = FALSE])
plot(sac_knn_m)
# Overlay analytical curves for comparison
Use the groups argument to run a single accumulation
with automatic native/alien splitting. Using 1990s subset.
rds_grouped <- file.path(comp_dir, "grouped_status.rds")
if (file.exists(rds_grouped)) {
sac_grouped <- readRDS(rds_grouped)
} else {
sac_grouped <- spacc(all_mat[idx_90s, , drop = FALSE], coords[idx_90s, ],
n_seeds = 50, method = "knn",
distance = "haversine", groups = status_vec, seed = 42)
saveRDS(sac_grouped, rds_grouped)
}
plot_grouped_normalized(sac_grouped,
title = "Grouped by Status (normalized)")
Compare accumulation curves across broad habitat types (EUNIS level 1). We subset to habitats with at least 100 plots.
habitat <- header$Eunis_lvl1
hab_counts <- table(habitat)
keep_habs <- names(hab_counts[hab_counts >= 100])
keep_habs <- keep_habs[!keep_habs %in% c("", "~")] # drop junk codes
cat(sprintf("Habitats with >= 100 plots: %s\n", paste(keep_habs, collapse = ", ")))
## Habitats with >= 100 plots: P, Q, R, S, T, U, V
cat("Plot counts:\n")
## Plot counts:
print(hab_counts[keep_habs])
## habitat
## P Q R S T U V
## 223 1102 14193 1109 7001 327 4642
rds_hab <- file.path(comp_dir, "habitat_curves.rds")
if (file.exists(rds_hab)) {
hab_results <- readRDS(rds_hab)
} else {
hab_results <- lapply(keep_habs, function(h) {
idx <- which(habitat == h)
spacc(all_mat[idx, , drop = FALSE],
coords[idx, ],
n_seeds = 50, method = "knn",
distance = "haversine", seed = 42)
})
names(hab_results) <- keep_habs
saveRDS(hab_results, rds_hab)
}
hab_combined <- do.call(c, hab_results)
plot(hab_combined) + ggtitle("Accumulation by EUNIS habitat (level 1)")
Are aliens more concentrated in certain habitats?
hab_alien <- header[, .(PlotObservationID, Eunis_lvl1)]
hab_alien <- merge(hab_alien,
species[, .(n_native = sum(STATUS == "native"),
n_alien = sum(STATUS == "neo")),
by = PlotObservationID],
by = "PlotObservationID")
hab_summary <- hab_alien[Eunis_lvl1 %in% keep_habs,
.(mean_native = mean(n_native),
mean_alien = mean(n_alien),
alien_pct = mean(n_alien / (n_native + n_alien)) * 100),
by = Eunis_lvl1]
print(hab_summary[order(-alien_pct)])
## Eunis_lvl1 mean_native mean_alien alien_pct
## <char> <num> <num> <num>
## 1: V 18.189573 2.65812150 14.7518372
## 2: U 15.547401 0.57186544 2.7899982
## 3: Q 11.684211 0.25045372 2.1949699
## 4: T 34.986288 0.65119269 1.9088159
## 5: R 29.732403 0.50806736 1.8607575
## 6: S 26.117223 0.23805230 0.9541968
## 7: P 5.834081 0.01793722 0.1002564
rds_hab_na <- file.path(comp_dir, "habitat_native_alien.rds")
if (file.exists(rds_hab_na)) {
hab_na_results <- readRDS(rds_hab_na)
} else {
hab_na_results <- lapply(keep_habs, function(h) {
nm <- native_mat[which(habitat == h), , drop = FALSE]
am <- alien_mat[which(habitat == h), , drop = FALSE]
nm <- nm[, colSums(nm) > 0, drop = FALSE]
am <- am[, colSums(am) > 0, drop = FALSE]
if (ncol(nm) < 2 || ncol(am) < 2 || nrow(nm) < 3) return(NULL)
idx <- which(habitat == h)
list(
native = spacc(nm, coords[idx, ],
n_seeds = 30, method = "knn",
distance = "haversine", seed = 42),
alien = spacc(am, coords[idx, ],
n_seeds = 30, method = "knn",
distance = "haversine", seed = 42)
)
})
names(hab_na_results) <- keep_habs
saveRDS(hab_na_results, rds_hab_na)
}
for (h in keep_habs) {
if (is.null(hab_na_results[[h]])) next
print(plot_normalized(Native = hab_na_results[[h]]$native,
Alien = hab_na_results[[h]]$alien,
title = sprintf("Habitat %s: Native vs Alien (normalized)", h)))
}
Track phylogenetic diversity (MPD, MNTD) as the sampling window expands spatially. Using 1990s subset for speed.
library(ape)
tree_cache <- file.path(data_path, "phylo_tree.rds")
if (file.exists(tree_cache)) {
cat("Loading cached phylogeny...\n")
tree <- readRDS(tree_cache)
run_phylo <- TRUE
} else if (!requireNamespace("V.PhyloMaker2", quietly = TRUE)) {
cat("V.PhyloMaker2 not installed and no cached tree found - skipping.\n",
"Install with: remotes::install_github('jinyizju/V.PhyloMaker2')\n")
run_phylo <- FALSE
} else {
run_phylo <- TRUE
library(V.PhyloMaker2)
sp_lookup <- unique(species[, .(species = WFO_TAXON, family = WFO_FAMILY)])
sp_lookup[, genus := sub(" .*", "", species)]
cat(sprintf("Placing %d species on mega-tree (first run, will cache)...\n",
nrow(sp_lookup)))
phylo_result <- phylo.maker(
sp.list = as.data.frame(sp_lookup[, .(species, genus, family)]),
scenarios = "S3"
)
tree <- phylo_result$scenario.3
tree$tip.label <- gsub("_", " ", tree$tip.label)
saveRDS(tree, tree_cache)
cat(sprintf("Tree saved to %s\n", tree_cache))
}
## Loading cached phylogeny...
if (run_phylo) {
# Reload matrices for subsetting (freed above)
mats <- readRDS(file.path(comp_dir, "matrices.rds"))
all_mat_phylo <- mats$all; native_mat_phylo <- mats$native; alien_mat_phylo <- mats$alien
rm(mats); invisible(gc())
in_tree <- colnames(all_mat_phylo) %in% tree$tip.label
cat(sprintf("Species matched to tree: %d / %d (%.0f%%)\n",
sum(in_tree), ncol(all_mat_phylo), mean(in_tree) * 100))
phylo_mat <- all_mat_phylo[idx_90s, in_tree, drop = FALSE]
tree <- keep.tip(tree, colnames(all_mat_phylo)[in_tree])
}
## Species matched to tree: 2697 / 2700 (100%)
if (run_phylo) {
rds_phylo <- file.path(comp_dir, "phylo_sac.rds")
if (file.exists(rds_phylo)) {
phylo_sac <- readRDS(rds_phylo)
} else {
phylo_sac <- spaccPhylo(phylo_mat, coords[idx_90s, ], tree = tree,
metric = c("mpd", "mntd"),
n_seeds = 30, method = "knn",
distance = "haversine", seed = 42)
saveRDS(phylo_sac, rds_phylo)
}
print(phylo_sac)
plot(phylo_sac)
}
## spacc phylogenetic diversity: 8814 sites, 2697 species, 30 seeds
## Metrics: mpd, mntd
if (run_phylo) {
rds_phylo_na <- file.path(comp_dir, "phylo_native_alien.rds")
if (file.exists(rds_phylo_na)) {
pna <- readRDS(rds_phylo_na)
phylo_native <- pna$native; phylo_alien <- pna$alien
rm(pna)
} else {
native_in_tree <- colnames(native_mat_phylo) %in% tree$tip.label
alien_in_tree <- colnames(alien_mat_phylo) %in% tree$tip.label
phylo_native <- spaccPhylo(native_mat_phylo[idx_90s, native_in_tree, drop = FALSE],
coords[idx_90s, ], tree = tree, metric = "mpd",
n_seeds = 30, method = "knn",
distance = "haversine", seed = 42)
phylo_alien <- spaccPhylo(alien_mat_phylo[idx_90s, alien_in_tree, drop = FALSE],
coords[idx_90s, ], tree = tree, metric = "mpd",
n_seeds = 30, method = "knn",
distance = "haversine", seed = 42)
saveRDS(list(native = phylo_native, alien = phylo_alien), rds_phylo_na)
}
cat("Native MPD accumulation:\n"); print(phylo_native)
cat("\nAlien MPD accumulation:\n"); print(phylo_alien)
}
## Native MPD accumulation:
## spacc phylogenetic diversity: 8814 sites, 2375 species, 30 seeds
## Metrics: mpd
##
## Alien MPD accumulation:
## spacc phylogenetic diversity: 8814 sites, 322 species, 30 seeds
## Metrics: mpd
Track functional diversity (FDis) as sites accumulate. Using family + status as proxy traits. Using 1990s subset.
if (!requireNamespace("FD", quietly = TRUE)) {
cat("FD package not installed - skipping functional analysis.\n",
"Install with: install.packages('FD')\n")
run_func <- FALSE
} else {
run_func <- TRUE
# Reload all_mat for subsetting
mats <- readRDS(file.path(comp_dir, "matrices.rds"))
all_mat_func <- mats$all
rm(mats); invisible(gc())
sp_info <- unique(species[, .(species = WFO_TAXON, family = WFO_FAMILY,
status = STATUS)])
sp_info <- sp_info[species %in% colnames(all_mat_func)]
sp_info <- sp_info[!duplicated(species)]
setkey(sp_info, species)
trait_df <- data.frame(
family = as.factor(sp_info$family),
status = as.factor(sp_info$status),
row.names = sp_info$species
)
shared_sp <- intersect(colnames(all_mat_func), rownames(trait_df))
func_mat <- all_mat_func[idx_90s, shared_sp, drop = FALSE]
trait_df <- trait_df[shared_sp, , drop = FALSE]
# Convert factors to numeric dummy variables for distance computation
trait_mat <- model.matrix(~ family + status - 1, data = trait_df)
rownames(trait_mat) <- rownames(trait_df)
cat(sprintf("Functional analysis: %d species with traits, %d sites\n",
length(shared_sp), nrow(func_mat)))
rm(all_mat_func); invisible(gc())
}
## Functional analysis: 2699 species with traits, 8814 sites
if (run_func) {
rds_func <- file.path(comp_dir, "func_sac.rds")
if (file.exists(rds_func)) {
func_sac <- readRDS(rds_func)
} else {
func_sac <- spaccFunc(func_mat, coords[idx_90s, ], traits = trait_mat,
metric = "fdis",
n_seeds = 30, method = "knn",
distance = "haversine", seed = 42)
saveRDS(func_sac, rds_func)
}
print(func_sac)
plot(func_sac)
}
## spacc functional diversity: 8814 sites, 2699 species, 124 traits, 30 seeds
## Metrics: fdis
Use the halo dataset (Austria + surrounding buffer) to see how edge effects influence accumulation.
# Load areaOfEffect countries data (needed for support = "Austria")
data(countries, package = "areaOfEffect", envir = environment())
header_halo <- fread(file.path(data_path, "austria_header_halo.csv"))
species_halo <- fread(file.path(data_path, "austria_species_halo.csv"))
cat(sprintf("Total plots: %d (core: %d, halo: %d)\n",
nrow(header_halo),
sum(header_halo$aoe_class == "core"),
sum(header_halo$aoe_class == "halo")))
## Total plots: 116413 (core: 52546, halo: 63867)
common_halo <- intersect(header_halo$PlotObservationID,
species_halo$PlotObservationID)
header_halo <- header_halo[PlotObservationID %in% common_halo]
coords_halo <- header_halo[, .(x = Longitude, y = Latitude)]
plot_ids_halo <- header_halo$PlotObservationID
halo_mat <- make_pa_matrix(species_halo[, .(PlotObservationID, WFO_TAXON)],
plot_ids_halo)
rds_halo <- file.path(comp_dir, "halo.rds")
if (file.exists(rds_halo)) {
h <- readRDS(rds_halo)
sac_halo <- h$halo; sac_core <- h$core
rm(h)
} else {
sac_halo <- spacc(halo_mat, coords_halo, n_seeds = 30, method = "knn",
distance = "haversine", support = "Austria",
include_halo = TRUE, seed = 42)
sac_core <- spacc(halo_mat, coords_halo, n_seeds = 30, method = "knn",
distance = "haversine", support = "Austria",
include_halo = FALSE, seed = 42)
saveRDS(list(halo = sac_halo, core = sac_core), rds_halo)
}
halo_combined <- c(`With halo` = sac_halo, `Core only` = sac_core)
plot(halo_combined)
sessionInfo()
## R version 4.5.2 (2025-10-31 ucrt)
## Platform: x86_64-w64-mingw32/x64
## Running under: Windows 11 x64 (build 26200)
##
## Matrix products: default
## LAPACK version 3.12.1
##
## locale:
## [1] LC_COLLATE=English_United States.utf8
## [2] LC_CTYPE=English_United States.utf8
## [3] LC_MONETARY=English_United States.utf8
## [4] LC_NUMERIC=C
## [5] LC_TIME=English_United States.utf8
##
## time zone: Europe/Luxembourg
## tzcode source: internal
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] ape_5.8-1 hexify_0.3.8 rnaturalearth_1.2.0
## [4] sf_1.0-24 ggplot2_4.0.1 data.table_1.18.0
## [7] spacc_0.1.0
##
## loaded via a namespace (and not attached):
## [1] magic_1.6-1 sass_0.4.10 generics_0.1.4
## [4] class_7.3-23 KernSmooth_2.23-26 lattice_0.22-7
## [7] digest_0.6.39 magrittr_2.0.4 evaluate_1.0.5
## [10] grid_4.5.2 RColorBrewer_1.1-3 fastmap_1.2.0
## [13] Matrix_1.7-4 jsonlite_2.0.0 e1071_1.7-17
## [16] DBI_1.2.3 gridExtra_2.3 mgcv_1.9-3
## [19] viridisLite_0.4.2 scales_1.4.0 permute_0.9-8
## [22] ade4_1.7-23 jquerylib_0.1.4 abind_1.4-8
## [25] cli_3.6.5 rlang_1.1.7 units_1.0-0
## [28] splines_4.5.2 withr_3.0.2 cachem_1.1.0
## [31] yaml_2.3.12 vegan_2.7-2 otel_0.2.0
## [34] geometry_0.5.2 tools_4.5.2 parallel_4.5.2
## [37] dplyr_1.1.4 vctrs_0.7.0 R6_2.6.1
## [40] proxy_0.4-29 lifecycle_1.0.5 classInt_0.4-11
## [43] MASS_7.3-65 cluster_2.1.8.1 pkgconfig_2.0.3
## [46] RcppParallel_5.1.11-1 pillar_1.11.1 bslib_0.9.0
## [49] gtable_0.3.6 glue_1.8.0 Rcpp_1.1.1
## [52] rnaturalearthdata_1.0.0 xfun_0.55 tibble_3.3.1
## [55] tidyselect_1.2.1 FD_1.0-12.3 knitr_1.51
## [58] farver_2.1.2 nlme_3.1-168 htmltools_0.5.9
## [61] rmarkdown_2.30 labeling_0.4.3 compiler_4.5.2
## [64] S7_0.2.1